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Abstract 

We use tools of the equilibrium statistical mechanics of disordered systems 
to study analytically the statistical properties of an ecosystem composed of 
N species interacting via random, Gaussian interactions of order p > 2, and 
deterministic self- interactions u > 0. We show that for nonzero u the effect 
of increasing the order of the interactions is to make the system more coop- 
erative, in the sense that the fraction of extinct species is greatly reduced. 
Furthermore, we find that for p > 2 there is a threshold value which gives 
a lower bound to the concentration of the surviving species, preventing then 
the existence of rare species and, consequently, increasing the robustness of 

the ecosystem to external perturbations. 
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Conservationists' arguments in favor of biodiversity have often appealed to the existence 
of intricate ties among apparently unrelated species in which, for instance, the strengths 
of the interactions between any pair of species would depend on the concentrations or fre- 
quencies of a variety of different ones [1]. Although the roles played by the total number 
of species, as well as by the strengths of their pairwise interactions, in the stability of an 
ecosystem are now fairly well understood both theoretically and experimentally H-§J, it is 
still not clear whether high-order interactions among the species would actually bring any 
advantage, in the sense of a larger robustness, to the ecosystem. This is the main issue 
we address in this letter. Albeit the model proposed here is somewhat unrealistic from 
the biological viewpoint, since its dynamics is governed by a Lyapunov function, it clearly 
points out the advantages of high-order interactions, making clear-cut, nontrivial predictions 
such as the reduction of the number of extinct species and the existence of a concentration 
threshold which excludes rare species from the ecosystem at equilibrium. 

Traditionally, the study of co-evolution of species has been restricted to deterministic 
interactions |5[], however the ever-present uncertainties about how the species are actually 
interacting allied to the overwhelming complexity of those interactions [Q motivate an al- 
ternative, and perhaps complementary, approach in which the strengths of the interactions 
between the species are assigned at random. In this letter we solve analytically a model 
of co-evolution of iV species interacting via random, high-order interactions. Our model is 
a generalization of the random replicant model |6|-[8| which considers only pair interactions 
between the species. 

We consider an infinite population (ecosystem) composed of individuals belonging to N 
different species whose fitness Ti [i = 1,...,N) are the derivatives Ti = dT/dxi of the 
fitness functional T defined as 

-F T~~lp (x) U ^ ] Xj + ^ t Ji\i2...i v %i\%i2 ■ ■ ■ %i p (1) 

i l<ii<i2-..<i p <N 

where Xi/N is the concentration of species i. These variables satisfy the constraints 
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and Xi > Vi. Here the coupling strengths are statistically independent random variables 
with a Gaussian distribution 
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for zi < Z2 < . . . < i p . The self-interaction parameter u > acts as a cooperation pressure 
limiting the growth of any single species, and it is crucial to guarantee the existence of a 
nontrivial thermodynamic limit, N — > oo. It can be shown that the dynamics 

dxj 

dxi N ^ ~ K dx k 



dt Xl 
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minimizes Ti p (x) while the mean Y^i%i is a constant of motion (see, e.g., Ref. page 240). 
This type of first-order differential equation, termed replicator equation, has been used to 
describe the evolution of self-reproducing entities (replicators) in a variety of fields, such as 
game theory, prebiotic evolution and sociobiology, to name only a few ||. In particular, 
a fourth-order interactions replicator equation was shown to govern the game dynamics in 
Mendelian (sexual) populations JTD|. 

For the sake of simplicity, in writing the fitness functional, Eq. ([!]), we have implicitly as- 
sumed that the couplings Ji 1 i 2 „.i J) are invariant under permutations of the indices h, . . . ,i p . 
We must stress, however, that regardless whether the couplings are invariant or not, the 
interaction term in the replicator equation, namely &H p /dxi, will be invariant to permu- 
tations of the species indices, and so the dynamics will converge to a fixed point. In this 
sense, the mere existence of a fitness functional (Lyapunov function) is a severe assumption 
from the biological viewpoint. On the other hand, it allows full use of the tools of the equi- 
librium statistical mechanics to study analytically the properties of the fixed points of the 
corresponding replicator equation. 

In the sequel we present the results of the replica analysis of the statistical properties 
of the ground state of the multispecies interaction Hamitonian (|l|). Following the standard 



prescription of performing quenched averages on extensive quantities only |TT], we define 
the average free-energy density / as 

-Pf= lim \-A\nZ) (5) 

where 

2 = J™ IT dx i 6 (n - £ x )j e"^W (6) 

is the partition function and ft — 1/T is the inverse temperature. Taking the limit T — ► 
in Eq. ([|) ensures that only the states that minimize 7i p (x) will contribute to Z. Here (. . .) 
stands for the average over the coupling strengths. As usual, the evaluation of the quenched 
average in Eq. (|5]) can be carried out through the replica method [11]. Within the replica- 



symmetric framework we find that, in the thermodynamic limit, the average ground-state 
energy per species is given by 



6 = lim / = u C Dzxl{z)- V -y q^ 1 (7) 
where Dz = dzexp (—z 2 /2) j\phx is the Gaussian measure, 
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and x s (z) is the positive solution of 
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-p (p - 1) yq p ~ 2 x s - puxl- 1 - \^f- x ) (A + z) = 0, (9) 
which maximizes the effective Hamiltonian 

E x = ^ P (p - 1) y<f ~ V - «x p - (^q p ~ l ) V2 (A + z) x. (10) 

We note that x s (z) = for z > 7. Here the saddle-point parameters g, and A are given 
by the equations 

1= r Dzx s (z), (11) 

J — 00 



q = f Dzx 2 s ( z ) 

J — oo 
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Although in general these equations can be solved numerically only, we can easily obtain an 
analytical solution for large u: 
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The physical order parameter q is defined by [[LI]] 
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where (• • -)t stands for a thermal average taken with the Gibbs probability distribution 



W (x) = I <5 [ JV - $>i ) exp (x) 



(18) 



Hence, values of g of order of 1 indicate the coexistence of a macroscopic number of species, 
while large values of q signal the dominance of a few species (i.e., the number of surviving 
species increases like N x with x < 1) only. In Fig. |l]we present the physical order parameter 
q as a function of the cooperation pressure u for several values of p. As expected, for large u 
the ecosystem is cooperative, in the sense that almost all species survive, and so q ~ 1. For 
small u the system enters a strongly competitive regime characterized by the divergence of q, 
though the onset of this regime can be postponed by increasing the order of the interactions 
p, as illustrated in the figure. Interestingly, the analysis of the effective Hamiltonian (|ToD 
for u = shows that x s and consequently q [see Eq. (0)] are finite only for p < 1, which 



corresponds to a random version of Szathmary's model of parabolic growth [12,13 1 . As 
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the divergence of q signals the survival of only a few species, the finitude of q at u = is 
consistent with a parabolic growth for which the coexistence of all species is assured. 

To better understand the distribution of species in the ground state we calculate the 
distribution of probability that a certain species concentration, say Xk, assumes the value x, 
defined by 



with W (x) given by Eq. (|T^). As all species concentrations are equivalent we can write 
Vk {x) = V [x] Vfc. Moreover, to handle a possible singularity in the limit (3 — > oo it is more 
convenient to consider instead the cumulative distribution function C (x) = Jq dx' V (x') . 
Carrying out the calculations we obtain 



where Q(x) = 1 for x > and otherwise. For u^oowe find C(x) — (x — 1) regardless 
of the value of p since in this case the equilibrium solution is = 1 Vi. An interesting 
feature of the cumulative distribution function is that C (0) is nonzero, indicating thus that 
the probability distribution V (x) has a delta peak at x = 0. In fact, the first term of the rhs 
of Eq. i.e., C (0), yields the fraction of extinct species in the ground state. Moreover, 
as shown in Fig. the constancy of C(x) up to a threshold concentration value x t = x s (7) 
indicates that there is a lower bound to the concentration of any surviving species. As 
illustrated in Fig. [|, x t decreases with the cooperation pressure u, and increases with the 
order p of the interactions. Since x t = at any finite value of u for p = 2, a nonzero 
value of Xt is an effect of the higher-order of the interactions. Clearly, the nature of this 
concentration threshold is totally distinct from that of the threshold obtained in the limit 
u — > 00, which equals 1 for all p. As expected, x t — > 00 for u = since C (x) — 1 for all x in 
this limit, while 
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for u large. Of course, the existence of such a threshold has far-reaching consequences on 
the stability and robustness of the population against external perturbations since it implies 
the absence of rare, and hence proner to extinction, species whose loss might cause dramatic 
effects in the whole population 0. Furthermore, for fixed u the fraction of extinct species 



which shows that increasing the order of the interactions among the species makes the 
ecosystem more cooperative, thus corroborating the conclusions drawn from the analysis of 



We have verified the validity of the replica-symmetric solution by performing the standard 
stability analysis ||14|| . In particular, that solution becomes unstable for u smaller than 
1/V2 « 0.707 and 0.106 for p = 2 and 3, respectively, while for p > 4 we find that it is 
unstable only at u — 0. Hence our main results are not affected by the (local) instability of 
the replica-symmetric solution. 

To conclude we must emphasize that our results describe the equilibrium properties of the 
population only. Important issues such as whether the absence of rare species at equilibrium 
would imply that the ecosystem is stable with respect to the invasion of rare mutant species, 
or whether the effect of a perturbation decreasing the concentration of a single species to 
a value below x t would lead to the collapse of the entire ecosystem, can be addressed only 
through a dynamical approach ||, which is beyond the scope of our present work. We hope 
the nontrivial predictions of our model will provide motivation for the proposal of more 
realistic models of high-order multispecies interactions. 

We thank Peter F. Stadler and Rita M. Z. dos Santos for useful discussions. The work 
of J.F.F. is supported in part by Conselho Nacional de Desenvolvimento Cientifico e Tec- 
nologico (CNPq) and Fundagao de Amparo a Pesquisa do Estado de Sao Paulo (FAPESP), 
Proj. No. 99/09644-9. V.M.O. is supported by FAPESP. 



decreases with increasing p (see Fig. [|). More pointedly, for large u we find 
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Fig. |. 
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FIGURES 

FIG. 1. Physical order parameter q as a function of the cooperation pressure u for p = 2, 3, 5 
and 10. 

FIG. 2. Cumulative distribution function of the ground-state species concentrations for p = 3 
and several values of u as indicated in the figure. The dashed line is the result for u — > oo. 

FIG. 3. Concentration threshold xt as a function of u for p = 3, 5, 7, 9, 11 and 13. Note that 
Xt = for p = 2. 

FIG. 4. Fraction of extinct species in the ground-state C(0) as a function of u for several values 
of p as indicated in the figure. Note that C(0) = 1 for u = indicating that only a few species 
survive in the thermodynamic limit. 
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